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ABSTRACT 

In this paper we present calculations of sound radiated from unflanged 
cylindrical ducts. The numerical simulation models the problem of an aero- 
engine inlet. The time-dependent linearized Euler equations are solved from a 
state of rest until a time harmonic solution is attained. A fourth-order 
accurate finite difference scheme is used. Solutions are obtained from a 

fully vectorized Cyber-203 computer program. Cases of both plane waves and 

spin modes are treated. Spin modes model the sound generated by a turbofan 
engine. Boundary conditions for both plane waves and spin modes are treated. 
Solutions obtained are compared with experiments conducted at NASA Langley 
Research Center. 
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IMTRODUCTIOH 


In this paper we present a con^utatlonal method to study sound radiated 
from an unflanged cylindrical duct. An incident field which is either a plane 
wave or a spinning mode (i.e., dependence on the azimuthal angle) propagates 
do\m the duct. At the open end of the duct, sound is radiated out into the 
farfleld and a reflected wave traveling upstream in the duct is generated. 
This problem is of importance in the study of noise radiated from aero-engine 
Inlets and in the development of effective duct liners. 

A significant amount of work has been done on the computation of sound 
propagation in an infinitely long duct. A survey of such work may be found in 
[1]. The open end of the duct and the ensuing outward radiation of energy 
significantly con^jllcates the problem. A further complication is the presence 
of the inlet flow about which little is known experimentally. In the 
procedure adapted here the solution Is obtained by solving the Euler equations 
linearized about an arbitrary mean flow. Thus the method is general enough to 
permit comiputation of the linearized fluctuating field about an experimentally 
determined mean flow. However, in this paper only the case of no mean flow is 
considered . 

We will briefly discuss some work which has been done in the past and is 
relevant to our work. The earliest work in calculating the sound wave 
(pressure) radiated from cylindrical ducts is due to Levine and Schwinger [8] . 
They provided a method to predict sound from a semi-infinite thin pipe, when a 
plane wave is incident upstream in the pipe, using the Welner-Hopf technique. 
This work motivated several other researchers in this field, in particular 
Savkar [10] provided a method to predict sound using Weiner-Hopf techniques 
for the case of an Incident spinning mode. Ting and Keller [13] developed an 
asymptotic expansion valid for plane wave incidence and flow frequencies. For 
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higher frequencies asymptotic methods have not been successfully applied to 
this problem because there are different length scales inside the pipe and in 
the farfleld. Though these methods provide some means to compute the sound 
radiated from engine Inlets, they do not correspond to the entire physical 
situation due to the thickness of the inlets. With a given thickness of the 
duct and for smooth geometries, calculations are effectively handled by 
integral equation methods. One such work is due to Horowitz, et al. [5], 
This method is based on a Helmholtz equation approach for the fluctuating 
sound pressure field and does not have the means to incorporate a mean flow. 
Thus a method is needed to study the radiation of sound that has the 
capability of handling flows. This is the motivation of this paper. 


2. FORMULATION OF THE PROBLEM 


We obtain the field equations from the fluid equations. The equations of 
Invlscid flow with the standard summation convention can be written as a first 
order system 

+ div(pv) = 0 



( 2 . 1 ) 


Here P is the density, v the velocity, and p is the pressure. We divide 
the flow variables into mean and fluctuating parts. We thus write 


P = P + p' 

V = U + u' (2.2) 

p = p + p' , 
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where the bar denotes a mean quantity Independent of time. We reformulate the 
resulting system by replacing the fluctuating density p' by the fluctuating 
pressure p' which Is the common acoustic variable. We assume that the flow 
is homentroplc and has no mean temperature gradient. It then follows that 

p = Ap^ (2.3) 

or 

p' ~ ^ 0(p (2.4) 

where cq Is the ambient speed of sound. Thus the resulting system from 
(2.1) becomes 


1 . 1 


2 ' 2 dlv(Pii'') -div(p U) + q 


(2.5) 


"(it + 


j 


3U. , 3U. 

+ u^ — ^ U — — 

T Uj J T ^ Ij. ^ 


j 3x 


J 


'0 


4- i£l 

j “"j s 


a- 3U. 

^ 3l^ ‘I’ 


where q denotes higher-order terms containing p p u , u etc. The 
left hand of system (2.5) contain the first order interacting terms between 
the fluctuating and mean quantities. The right hand side contain the mean 
flov7 and fluctuating quantities of the lower order terms. In general the 
system (2.5) is a linear first-order hyperbolic system which includes all of 
the first-order terms for the fluctuating field subject to a time-dependent 
inflow condition. 

For the present purpose we assume the mean flow is zero. Thus the system 


(2.5) reduces to 
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4 PQ ° 

( 2 . 6 ) 

3u' 

PQ 

where Pq Is the density of the ambient fluid. We non-dimensionalize these 

equations. Length is non-dimensionalized by the diameter of the pipe (d), 

2 

time by Cg/d, pressure by PqCq and the velocity by Cq to obtain 


i£ 

3t 


+ div u = 0 


3u 

^ + Vp = 0. 


(2.7) 


Note that in (2.7) the prime on the fluctuating quantities are dropped for 
convenience. 

REMARK 2.1 ; If p and u are time harmonic, that is 

p(x,t) = p(x)e“^*^’^. 


u(x,t) = u(x)e~^^^ 


where k is the wave number then the system (2.7) reduces to 


Ap + k^p = 0. (2.8) 

The problem discussed here is to solve the system (2.7) for p and u 
subject to appropriate boundary conditions which will be discussed later. 
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The technique here is to drive the system (2.7) with a time harmonic 
source which will yield the time harmonic solutions, namely the solution of 
(2«8). This technique is essentially the numerical implementation of an 
appropriate limiting amplitude principle. For exterior problems this method 
has been demonstrated by Krlegsmann and Morawetz [6] and Taflove and 
Umashankar [12] and for wave guide problems by Baumeister [1] and Krlegsmann 
[7]. 


3. SOLUTION ^OCEDURES 


We take the origin at the open end of the duct so that its generators are 
parallel to the z axis (Figure 1). We look for solutions of (2.7) of the 
form 

p(r,0,z,t) = P(r,z,t)e^”^ 


u(r,0,z,t) 


U(r,z,t)e 


im0 


(3.1) 


In general solution will have the form 


= I 

m“0 


P^(r,z,t)e 


im0 


but we confine our solutions as In (3.1) for a single mode m. For m > 1 
the solution is referred to as the spin mode solution. In physical situations 
they correspond to the modes provided by a turbofan engine. For m = 0 (3.1) 

describes a plane wave. The more Important case is when m > 0 (a spinning 
mode). Then the system (2.7) becomes 
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3t 


+ u 

z 


+ V 

r 


+ 


V + Imw 
r 


= 0 


3u ^ _ 

3t 3z " 


0 


3v , 
at 3r 


= 0 


3w , 

3t r 


P 


= 0 


(3.2) 


and the problem reduces to solving the system (3.2) with appropriate boundary 
conditions. 

The method we use to solve (3.2) is an explicit method which is fourth- 
order accurate in space and second-order accurate in time and is due to 
Gottlieb and Turkel [4]. This is a higher-order accurate version of the 
MacCormack scheme. The use of an explicit method has been recently advocated 
by Baumeister [1]. The major advantages are drastically reduced storage 
requirements and programming ■ simplicity. Baumeister demonstrated the 
effectiveness of this approach for internal sound propagation with spinning 
modes [2]. The work in this paper extends this idea to the full radiation 
problem with a more accurate computational scheme. 

A typical computational domain is depicted in Figure 2. Referring to 
this figure, the computations are carried out in the rectangular region which 
is bounded by an inflow boundary and the farfield boundary. Thickness of the 
duct is allowed by a mesh size thickness in the r direction. The solution 
for large times is extremely sensitive to the inflow condition and also the 
farfield condition [9]. The solution is assumed to start from a state of 
rest, l.e., P = u = v = w = 0 at time t = 0. The system (3.2) can be 


written in the form 
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+ F + G = H , 
— t —2 — r — ’ 


(3.3) 


where 


w = 


P 


u 


V 


M 

_ V + imw 
r 

u 

, F == 

p 

0 

II 

0 

, and H = 

0 

V 


0 


P 


0 

w 


0 


0 


_ imP 



L. -J 


^ — 


r 


. (3.4) 


We use the method of time splitting to advance the solution from time t to 
t + 2At. If L^(At) and L^(At) denote symbolic solution operators to the 


one-dimensional equations 


w^ + F = H, 
t 2 1 


Wt + Gr = H2, 


(3.5) 


then the solution to (3.3) is advanced by the formula 


w(t + 2At) = L^(At) L^(At) L^(At) L^(At) w(t). 


(3.6) 


This procedure is second-order accurate in time. The fourth-order accuracy in 
space depends on formulating the difference scheme. For a one-dimensional 
system, i.e., for (3.5), we have 


Wi(t + At) = w^(t) (7F^ - 8F^^^ + F^^ 2 ) + 


w^(t + At) = j [w^(t) + w^(t + At) +-g^ [-7F^ + 8F^_^ - F^_ 2 ) + AtH^] , 


(3.7) 
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where denote F evaluated at w^, etc. This formula contains a forward 
predictor and a backward corrector. This is second-order accurate in space. 
One can formulate another variant which contains a backward predictor and a 
forward corrector which is also second-order accurate in space. In order to 
achieve fourth-order accuracy we alternate (3.7) and its variant in each time 
step. If there are N intervals with nodes at Zj^(i=0,l ,• • • ,N) then the 
predictor in (3.7) cannot be used at i = N-1 and at i = N and the 
corrector cannot be used at 1 = 0 and 1. Similar situations occur for the 
other variant too. At these points we extrapolate the fluxes using third- 
order extrapolations. For the right boundary we use the extrapolation 
formulae 


N+1 


= 4F — 6F 
N N-1 


+ 4F 


N-2 


- F 


N-3 


(3.8) 


N+2 


= 4F 


N+1 




4F 


N-1 


- F 


N-2’ 


and for the left boundary 


Fq = 4F^ - 6F2 + 4F3 - F^ 


F 


1 



6Fj^ + 4F2 - F3. 


(3.9) 


Since we are interested in time harmonic solutions, the numerical 
solution is monitored until the transient has passed out of the computational 
domain and the solution achieves a steady time harmonic dependence. 
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4. BOUNDARY CONDITIONS 

A very important feature of our work Is In obtaining appropriate boundary 
conditions. We derive our boundary conditions appropriate to an experiment 
carried out at NASA Langley Research Center [14] . The boundary conditions 
consists of two major parts. The first part is derivation of an Inflow 
condition which will model correctly the sound source. Next, we need accurate 
farfield boundary conditions which will simulate outgoing radiation. 


Inflow Conditions 

To derive Inflow boundary conditions we consider the time harmonic case 
and in particular equation (2.8). We look for spinning mode solutions of the 
form 


p = P(r,z)e 


lm0 


This yields 


j, 

r 


3 

3r 


3r^ 


(k" - 


0 , 

9z 


(4.1) 


If we separate variables by setting 

P(r,z) = f(r) g(z), 

we obtain 

f(r) = J (8r) and g(z) = 
m 


where 
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(4.2) 


and I Is to be determined. 

If a is the radius of the duct (here a = V 2 due to the non- 
diraensionalization) then the usual boundary condition on the pipe is the hard 
wall condition 


A 

8n 


0 


on r = a. 


This gives 

4— (j (3r)l =0 on r = a 

dr ^ m ■' 

or 

3a = X (n=0,l,2,**» ,) 

mn 


where X 's are the zeros of the functions J'(z). From (4.2), using the 
mn m / \ f t e> 

appropriate subscript corresponding to X^^ we have 


2 2 2 
(Jl ay = (ka)^ - X^ . 
mn mn 


(4.3) 


Definition: If ka > X , we Say the mode (m,n) is cut-on. 

Otherwise it is said to be cut-off. 

We now consider only cut-on modes. Then the solution of (4.1) has the 

form 


P(r,z) = I 

n=0 


±iz 


;/(ka)^ - 


mn 


J (-22- r) 

m'- a J 


a e 
n 


(4.4) 
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It is necessary to consider the case of a single cut“on mode propagating down 
the pipe. In this situation n = 0. Dropping the corresponding zero 
subscripts In (4.3) we consider the values of H given by 



Then the general solution Inside the pipe can be written as a combination of 
an incoming wave and a reflected wave. That Is 


p(r,0,z) = (e + R(k)e 


) 


m'" m a 
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(4.5) 


where R Is the reflection coefficient and Is also a function of the wave 

number k. Recalling that p(r,9,z,t) = p(r,9,z)e”^^*" and 

iiTi0 

p(r,0,z,t) = P(r,z,t)e^ we have 


P(r,z,t) 


iZ z -i£ z 

= (e + R(k)e ) 


J (X 
m^ m a' 


(4.6) 


The reflection coefficient R Is unknown. Thus we must eliminate R in 
(4«6). This Is accomplished by taking the time derivative of (4.6) and 
subtracting the times the z derivative to get 


but 


P 

t 



-2ik 


i.Z z 
m 
e 



r-i -Ikt 


P = - u from (3.2). 

z t 


Thus the above equation becomes 
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1 iZ z . 

(P +~u) = -2ik e “ J (X 

)!> t m tn a 

m 

We impose this boundary condition on the inflow boundary z = -L. We 
boundary condition on v at the Inflow by using (3.2), 


iz + iZ 

3t 8r 


as 


0 , 


together with 


to give 



Note that the coefficient of P here contains a singular term at 

However P contains a term (from (4.6)) of the form 

“^m^^m when r = 0 (II) is simply replaced by v = 0. 


Conditions on the Wall 

3 P 3P 

On the duct wall = 0, but 

8n 3r ’ 

3v . 3P n c /o o\ 

Tt 3r “ ^ from (3.2). 


(I) 

obtain 


(II) 

r » 0. 


(Ill) 


This implies v = 0 on the wall. We note that a general Impedance condition 
simulating an acoustic liner can be handled without difficulty. 
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Conditions on the Axis 

When m = 0 the system (3.2) has only three equations for p, u, and v. 
The first equation of (3.2) contains a term. Thus the boundary condition 

on the axis in this case is 

V = 0 on r = 0 (m = 0). (IV) 


When m = 1 the last equation of (3.2) gives 


+ = 0 . 

9t r 


Here P contains a term like J, ~) for z close to -L. Thus Is 

1, a dt 

nonzero. But from the first equation of (3.2) we have 


v+lw=0 on r=0 (m=l). 


(V) 


For m > 2 the first and the last equations of (3.2) give 

V = 0, w = 0 on r = 0 (m > 2). (VI) 

Farfleld Conditions 

Radiation conditions are applied at the farfleld boundaries. The 
development follows that given in [3]. Let R be the distance [r = /r + z ) 
from the origin to a point In the farfleld (see Figure 2). The condition we 
Impose here is the first member of a family of nonreflecting boundary 
conditions which are accurate as R This condition is 
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where 


IP + IP + P = 0 

3t 3R R ’ 


iz 

3R 


3P ^ 3P ^ 

-r~ cos ot + — sin a, 
3z 3r 


where a is the angle from the z axis to the farfleld point. Using the 
second and third equations of (3.2) we have 


3P 

3R 


Uj. cos a - sin a. 


Thus the radiation condition becomes 

-- (u cos a + V sin a) + -^ = 0. (VII) 

ot t R 

The conditions I through VII were used to obtain the results discussed in the 
next section. 

5. NUMERICAL RESULTS 

We computed the solutions with the details given in Sections 3 and 4 on a 
CDC Corp. Cyber-203 machine. The algorithm described above is almost totally 
vectorizable. For a very low frequency plane wave case the typical number of 
grid points in the r,z plane were 80 x 100. The inflow boundary was kept 
at z = -lOd (10 diameters) and the radiation boundary was chosen so as to 
enclose a circle of radius lOd. For high frequencies the typical grid sizes 
were 115 x 135 and the Inflow boundary as varied from z = -lOd to z = - 8 d. 

To verify the effectiveness of the code we compared our results with 
asymptotic expansion obtained by Ting and Keller [13] for a low frequency 
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plane wave. To make comparisons we computed the solutions in the duct and on 
the axis at various stations for a non-dimensional frequency ka = 0.2. Here 

27f 

k = — is the wave number. Results are presented in Table I. 

For high frequencies we compared our results with the Weiner-Hopf results 
of Savkar and Edelfelt [11] and the experiment done at NASA Langley Research 
Center [14] . In this experiment the directivity patterns were measured on a 
circle at 10 diameters from the open end of the pipe. The test facility has a 
spin mode synthesizer which can produce both plane and spinning mode wave. 

The first comparison was made for the plane wave case (m = 0) and a 
non-dimensional frequency of ka = 3.76. Since the experimental results were 
obtained only in the farfleld the sound pressure level was plotted as a 
function of angle measured from z axis. The results are presented in Figure 
3 and show good agreement with the experiment. 

Figures 4 and 5 show typical comparisons of the spinning mode case 
with m = 2 and for frequency values ka = 3.37 and 4.40. As in these 
figures, except the plane wave case, the computed results agree within 5 dB 
levels. Clearly our results show better comparison than Welner-Hopf results 
[11] due to allowance of thickness. In these cases the results near the axis 
do not compare very well. This is due to the fact that in the experiment it 
is difficult to completely control other modes and plane waves. This is 
particularly true for this frequency ka = 4.40 which is close to the next cut- 
on mode. In the plane wave case the results were unexpectedly good. 
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6. VARIABLE GEOMETRY DUCTS 

We consider ducts with a local variable geometry cross section. The duct 
is assumed to be straight as z ->■ -«> (see Figure 4). Thus the Inflow 
boundary conditions previously formulated are still valid. The variable duct 
is incorporated in the numerical scheme by mapping it into a straight duct. 
This slightly changes the coefficients in the final system (3.3) but does not 
degrade the convergence to the time harmonic solution. 

Suppose the duct configuration is as in Figure 4. It has a curved 
boundary near z = 0 and has straight extension everywhere else. This allows 
us to have the same Inflow boundary conditions and the conditions on axis and 
also the radiation condition. But the boundary conditions on the wall will be 
changed. 

The Euler equations have the form 


w. +f +g +h=0. 
t z °r 


We do the following maps: 


( 6 . 1 ) 


( 6 . 2 ) 

r 

^1 an(z) 

We use chain rule to compute f^, in terms of f^^ and ate. where 

r = n(z) is the geometry of the duct. This yields 


w. + f 
t z. 


an(z) 


rn (z) ,'1 
2 ''r 

an(z) 1 


+ f 


n'(z) 

n(z) 


+ h = 0. 


(6.3) 
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This has the same form as (6.1). Thus very minor changes in the difference 
scheme and in the radiation boundary conditions are required. The boundary 

condition v = 0 on r = a Is replaced by the vanishing of the normal 

velocity on the wall. On the surface of the pipe a normal vector is (1, 

- an'(z)) in (r,z) coordinates. Thus the above condition reduces to 

V - an'(z)u = 0. (6.4) 

We simulated a duct where n(z) has the form 

( .5 z < -1 

n(z) = J 

( .5 - e(2z -l)(z+l) -1 < z < 0 

(see Figure 4) . The grids of the computational domain follow the same 
geometry. For e = .15 the results we obtained are shown in Figure 6 
compared with the straight duct situation. The dB level reduces at 90° by 

about 10 dB. This indicates the importance of the nozzle geometry in 

determining the farfield radiation pattern. 
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Table I. Comparison with Ting and Keller Solution 

ka = 0.2 


z 

Ting & Keller 
|P| 

Numerical 

|P| 

-10 

1.5026 

1.5054 

- 9 

1.0984 

1.1113 

- 8 

.3873 

.3544 

- 7 

.4355 

.3933 

- 6 

1.1133 

1.0874 

- 5 

1.7603 

1.7196 

- 4 

1.9280 

1.9583 

- 3 

1.8933 

1.9097 

- 2 

1.5495 

1.5477 

- 1 

1.0279 

1.0076 
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Figure 2 











dB, level 
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Figure 4 
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